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Abstract 

Background: Bone tissue is the main element of the human skeleton and is a 
dynamic tissue that is continuously renewed by bone-resorbing osteoclasts and 
bone-forming osteoblasts. 

The bone is also capable of repairing itself and adapting its structure to changes in 
its load environment through the process of bone remodeling. 
Therefore, this phenomenon has been gaining increasing interest in the last years 
and many laws have been developed in order to simulate this process. 

Results: In this paper, we develop a new law of bone remodeling in the context of 
damaged elastic by applying the thermodynamic approach in the case of small 
perturbations. 

The model is solved numerically by a finite difference method in the 
one-dimensional bone structure of a n-unit elements model. 

Conclusion: In addition, several numerical simulations are presented that confirm 
the accuracy and effectiveness of the model. 

Keywords: Numerical simulation, Thermodynamic approach, Small perturbations 
hypothesis, Bone remodeling, Bone density, Damage, Fatigue, Osteocyte, Elasticity, 
N-unit elements 



Introduction 

Bone is a living material that constantly replaces old tissue with new in a process called 
remodeling. It is also able to respond adaptively to its environment [1,2]. 

The bone remodeling process replaces approximately 20% of bone tissue annually; in 
healthy adults, bone remodeling occurs in a balanced, highly regulated manner in five 
phases: activation, resorption, reversal, formation, and quiescence as shown in Figure 1 
[3-5]. 

This process is assumed to repair the microdamage and maintain bone quality; and 
also occurs continuously with each cycle lasting 4 to 7 months [6]. 

Over the past, the progress made in understanding bone remodeling, through two 
groups: phenomenological approach and thermodynamic approach, has been truly 
spectacular. 

The thermodynamic approach was initiated first by chemists and was applied to con- 
tinuum mechanics by Eckart and Biot around 1950. Furthermore, this approach was 
adopted by introducing state variables [7,8] and thermodynamic potential which allows 
to define associated variables chosen for the study of the phenomenon [9]. 
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Figure 1 Bone remodeling sequence. 
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Many theoretical and experimental works have been published emphasizing the im- 
portance of applying this thermodynamical approach to a bone structure [10-13]. 

In this study, we adopt thermodynamic approach of irreversible processes to get a 
new law of bone remodeling taking into account the bone density described by the law 
of Mullender et al. [14] and the damage evolution in the bone proposed by Martin [15]. 

The found equation is solved by the finite difference method (FDM) in the one- 
dimensional bone structure of a n-unit elements model. 

Finally, we present some examples of numerical simulation results. 



Fundamental assumptions 

In order to construct a general framework for the description of the bone remodeling 
process, the following simplified assumptions are made: 

1) The bone is considered as a linear-elastic, isotropic and inhomogeneous material. 

2) The external mechanical loading acts as a stimulus for bone remodeling. 

3) The n-unit elements model is applied to the damaged-bone structure in the one- 
dimensional case. 

4) The small perturbations hypothesis (displacements and their time and spatial 
variations are small). 

5) The state coupling of damage with elastic strain. 

6) The remodeling processes can be considered isothermal, adiabatic and without 
internal generation of heat. 

7) The decoupling of the thermodynamic potential is assumed, such that: 
With: 

i// e {e e , D) is the thermodynamic potential depending on the elastic strain tensor e e and 
the damage variable D. 
i// r {(f>) is the thermodynamic potential depending on the bone density 0. 
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Thermodynamic approach 

The thermodynamics of irreversible processes allows the modeling of different materi- 
als behavior. This is accomplished by defining the state variables and the state potential 
and also the dissipation potential [16]. 

The general theory of adaptive damaged- elastic materials and general framework of 
continuum thermodynamics is considered to find a new law of bone remodeling 
[10,17]. 

A general definition of the thermodynamic forces associated with the internal vari- 
ables is given by: 

Where conjugate forces associated with internal variables (V^, k=l,2,. . .) by speci- 
fication of the thermodynamic potential if/(. . ., V^) as shown in Table 1 [17,18]. 

Within the hypothesis of small strains and small displacements, the state variables, 
observable and internal, are chosen in accordance with the physical mechanisms of 
deformation and degradation of the bone as follows [18]. 

Observable variables: 

- e is the total strain tensor associated with the stress tensor a. 

- T is the temperature associated with the specific entropy s. 

Internal variables: 

- e e is the elastic strain tensor associated with the stress tensor a. 

- D is the damage associated with a variable Y. 

- 0 is the bone density associated with the bone remodeling variable R. 

Table 2 summarizes the set of variables introduced [18]. 

We assume a bone remodeling variable, which is characterized by: 

R > 0 in the case of the formation phase. 
R = 0 in the case of the equilibrium phase. 
R < 0 in the case of the resorption phase. 

We postulate the existence of a thermodynamic potential from which the state laws 
can be derived [16,18,19]. 
The state potential: if/ = if/(e e , D, (p) 



Table 1 Thermodynamic variables 



State variables 




Associated 
variables 


Observable 


Internal 


T 




5 


8 




O 
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Table 2 Chart of thermodynamic variables 

State variables Associated 

Observable Internal variables 

T s 

£ O 

s e o 

D Y 

0 -R 



We assume the following decoupling: 

l^=^(8 e ,D) + ^(0) 

With: 

if/ e (e e ,D) is the thermodynamic potential depending on the elastic strain tensor and 
the damage variable. 
ij/ r {(j)) is the thermodynamic potential depending on the bone density. 
The associated variables are defined by: 

The second law of thermodynamics imposes a restriction on dissipation which can be 
represented in terms of the Clausius-Duhem inequality [7,16,18] : 

a:s-p.ij/>0 (3) 



The development of the Clausius-Duhem inequality: 

We note that: 

and 

(the thermo-elasticity law) [17] 
Inequality (3) may be written : 

Using the previous equations (1) (2) (4), we obtain: 
-Y.D+£.0>O (5) 



Discussion: 

In the case of constant damage: D = 0 
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The inequality (5) gives +R.0>O 

In the resorption area: 
0<O andRO 

=> +R.0>O 

In the formation area: 
0>O andR>0 

=> +R.0>O 

Equilibrium area (dead zone): 
0 = Oand + R.0 = O 

=> +R.0>O 

We can conclude that the inequality (5) is verified in accordance with the law of bone 
remodeling as shown in Figure 2 [20]. 

The conservation of energy equation: 

The first principle of thermodynamics gives: 

p.e = a : e e + r — divq (6) 



[17,21] 
With: 



- e is the specific internal energy. 

- r is the internal heat source. 

- q is the heat flux. 



We replaced pe by the expression derived from e = if/+ T.s [17,21] 



p.e = p.if/-\- p.T.s + p.T.s 



(7) 



d<t>/dt 





formation 



stimulus 



/ resorption 

Figure 2 Law of bone remodeling. 
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and p if/ by its expression as a function of the state variables and the associated 
variables. 

p.if/=a:e e +Y.D-R.(j) (8) 

We introduce the specific heat (capacity) defined by: 

C = T|| [16,17] and taking into account Fourier's law q = —k. gradT [17] 
We obtain: 

div q = —k.div(gradT) = —k.AT (9) 
Using the previous equations (6) (7) (8) (9), we can get: 
Y.D-R.(p + pT.s+p.T.s = r + kAT (10) 



Assumptions: 

The classical heat equation corresponds to a process: [17] 

- Without internal generation of heat created by the external sources: r=0. 

- With adiabatic evolution: k. AT= 0. 

- With isothermal transformation: T = 0 Therefore s = j . T = 0 Then s =0. 

Equation (10) may be written: 

Y.D-R.(p = 0 
Then, 

- D 

R = Y.- k (11) 

We have: Y = —Y 
Where 

Y = ^.E.e e :e e (12) 
(the strain energy release rate) [17] 

We have also: p.if/ e = ^.(1 — D).E.s e : s e (the strain energy) [17] and the equivalent 
constraint o eq is written by o eq = p^f = (1 — D).E.s e (The thermo-elasticity law) 
Then, 

(13) 



(1-D).E 
a eq = a in 1 

Using the previous equations (12) (13), equation (11) may be written: 



With: a eq = a in the one-dimensional case. 
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The Young's modulus of the bone which is an isotropic material and inhomogeneous 
is expressed as: 
E = (1 - D).E 0 with E Q = c.cff [22,23] 

c=100 and a =3 are two constants characteristic of the bone 
Then, 

E = (l-D).c.(p a 

Equation (14) may be written: 

a 2 D 1 

R((p,D) = . , . r 

W ; 2c (1-D) 3 (P a (P 

This equation represents the new law of bone remodeling developed by applying the 
thermodynamic approach in the context of damaged elastic. 

In this study, we introduce the law of damage developed by Martin [15] which shows 
that the damage in human cortical bone can grow exponentially until the fatigue failure 
[15,24,25]. 

The evolution law for the damage is expressed as: 

D = D 0 J dt Then, § = D =f d D 

With: 

f d : the fatigue life of the bone devoid of the remodeling [26] 

D 0 : the initial damage 

t : the time 

Finally, the new law of bone remodeling may be written as: 

J?(0 , D) = _^._O .J^ (15) 
W ^ 2c (1-D) 3 (P a (P V ^ 



Numerical resolution 

The new law of bone remodeling (Equation 15) was solved numerically by dividing it 
into three parts: 

1. a constant ^ 

2. a damage function 

3. a function of bone density ^ 

In this study, we use the law of bone density proposed by Mullender et al. [14]: 

90,. ^ <*k)(s k \ 

* r g^Ur v ) (16) 

With : 

— ^min — 4* — 4* max 

- (j) min is the density of completely resorbed bone 



Id ham mad and Abdali Theoretical Biology and Medical Modelling 2012, 9:51 
http://www.tbiomed.eom/content/9/1/51 



Page 8 of 1 1 



I* 



i" r 


r 


r 


r 




■ F/n . F " 1 


F„ 

I 


1 


? 

■ 


• 

3 

I 2 


• 

3 


t 

4 






nf 

m 


n 



I d 

N M 

total number of osteocytes 

Figure 3 Model geometry with n unit elements. 

I J 

- fimax is the maximum density defined for a compact bone 

- t is a positive constant related to the reaction time of bone tissue (constant of bone 
remodeling) 

- 1< i <n 

- (j)i density of bone tissue of element i 

- m (m < n) is the total number of osteocytes in the solid 

- Ik (1 < k <m) corresponds to the series of numbers of the elements containing an 
osteocyte 

- S k represents the density of deformation energy in I k 

- S re f reference stimulus value 

- is a parameter reflecting the intensity of the stimulus cell 

- d is the normalization factor limiting the area of influence of osteocyte 

- d(i,Ik) is the distance between the centers of geometric element i and the element 4 

We discretize into n-unit elements a bone fragment and we apply a compressive 
force evenly distributed over the various units (Figure 3) [23,27]. 

To solve the previous equation (16), we use the method of finite difference with an 
implicit scheme and the fixed point method [23,27,28]. 




Table 3 Values of the parameters used during the numerical simulations 



Data 


Symbol 


Values 


Units 


Maximum density 


0max 


1.75 


g/cm 3 


Minimal density 


0min 


0.01 


g/cm 3 


Initial density 


00 


0.6 


g/cm 3 


The step of time 


At 


5.1 0" 3 


UT 


The total force 


F 


10 


N 


The distance between 2 centers 


d 


25 


mm 


Reference stimulus value 


Sref 


0.04 


MPA 


The fatigue life of the bone 


f d 


3 


years 


n-unit elements of the bone fragment 




50 





m (m < n ) is the total number of osteocytes in the bone fragment. 
Constants a = 3 (3 = 0.5 D 0 = 0.8 c = 1 00 t = 1 . 
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We obtain: 

0° = 0O 



Results and discussion 

We simulated the case of a uniform distribution of the osteocyte cells, and of another 
heterogeneous case. The values of the parameters used during the numerical simula- 
tions are given in Table 3 [23,28]. 
Numerical results are shown in Figure 4. 



Evolution of the variable of bone remodeling 




c 

O 



C 

o 

-Q 



- simulation 6 

- simulation 5 

- simulation 4 
-simulation 3 
-simulation 2 

- simulation 1 



Evolution of the variable of bone remodeling 




Figure 4 Evolution of the variable of bone remodeling, (a) case of a uniform distribution of the 
osteocyte cells, (b) a heterogeneous case. 
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Figure 4 shows the temporal evolution of the variable of bone remodeling in the case 
of a uniform distribution of the osteocyte cells (n=m=50), and of another heteroge- 
neous case (n*m with m=30 in the central package of the osteocytes). 

The curves consists of three key periods. The first period of the curves corresponds 
to the resorption phase, where the variable of the bone remodeling was negative. The 
second period exhibits the formation phase, where the variable of the bone remodeling 
was positive. The third period defined as the interval between the resorption phase and 
the formation phase, which the variable of the bone remodeling reached zero. 

The resorption phase takes approximately 18 days, which is then followed by an equi- 
librium phase that can last for up to 10 days and finally by the formation phase from 
17 to 35 days. This is in agreement with results from the literature [6,29,30]. 

By comparing the curves in Figure 4 to the graph proposed by Terrier et al. [20], we 
see a good agreement. Furthermore, the curves that are found have a nonlinear shape. 

Conclusion 

In this paper, we proposed a thermodynamic approach in small perturbations for bone 
remodeling process. 

The adopted model takes into consideration both the bone density and the damage 
and gives a new law of bone remodeling. Then, the governing equation of the process 
was solved by the finite difference method in the one-dimensional bone structure with 
n-unit elements model. 

The numerical results obtained are in accordance with the experimental results found 
in the literature. 
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